Libraries and paths
library(easypackages)
libraries("here","tidyverse","ggeasy","patchwork","lmerTest","emmeans","reshape2")
plotdir = here("plots")
fontSize = 16
baseline_end = 15
transition_end = 35
# cols2use = c("#619cff","#156cff")
cols2use = c("#619cff","#ffa500")
Functions to use later
# functions to use for plotting
make_timeplot <- function(data4plot, x_var, y_var, grp_var,
baseline_end=15, transition_end=35,
yLabel,title2use,fontSize=16){
p = ggplot(data = data4plot,
aes_string(x = x_var, y = y_var, group = grp_var)) +
facet_grid(. ~ grp)
p = p + geom_smooth(se = FALSE,
colour='gray75',
alpha=0.1)
p = p + geom_smooth(se=TRUE, aes(group = interaction(grp,condition),
colour = condition,
fill = condition),
size=3)
p = p + geom_vline(xintercept = (baseline_end)) +
geom_vline(xintercept = transition_end) +
geom_hline(yintercept = 0) +
xlab("Time Window") +
ylab(sprintf("Baseline Normalized %s",yLabel)) +
guides(colour="none", fill="none") +
ggtitle(title2use) + easy_center_title() +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, size=fontSize),
strip.text.y = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize-2),
axis.title.x = element_text(size=fontSize),
axis.title.y = element_text(size=fontSize),
strip.text = element_text(size = fontSize))
return(p)
} # make_timeplot
# make_scatterplot
make_scatterplot <- function(data4plot, x_var, y_var, grp_var,
cols2use,xLabel,yLabel,title2use, dotSize=3){
p = ggplot(data = data4plot, aes_string(x = x_var, y = y_var, colour=grp_var)) +
geom_point(size = dotSize) +
geom_smooth(method=lm, se=FALSE) +
scale_colour_manual(values=cols2use) +
xlab(sprintf("Baseline Normalized %s",xLabel)) +
ylab(sprintf("Baseline Normalized %s",yLabel)) +
ggtitle(title2use) + easy_center_title() +
guides(colour="none") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, size=fontSize),
strip.text.y = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize-2),
axis.title.x = element_text(size=fontSize),
axis.title.y = element_text(size=fontSize),
strip.text = element_text(size = fontSize))
return(p)
} # function make_scatterplot
# make_blankplot
make_blankplot <- function(exp, dtype="H", lineSize = 3,ylimits =c(-6,6),fontSize = 16){
data_FR = data.frame(time = c(0,15,16,35,36,60, 0,15,16,35,36,60),
condition = c("baseline","baseline",
"transition","transition",
"treatment","treatment",
"baseline","baseline",
"transition","transition",
"treatment","treatment"),
grp = c("DREADD","DREADD","DREADD","DREADD","DREADD","DREADD",
"SHAM","SHAM","SHAM","SHAM","SHAM","SHAM"))
data_H = data_FR
if (dtype=="H"){
if (exp=="camk"){
data_FR$score = c(0,0,0,3,3,3, 0,0,0,0,0,0)
data_H$score = c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
title2use = "CamkII-hM3D(Gq)"
}else if (exp=="hm4di"){
data_FR$score = c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
data_H$score = c(0,0,0,3,3,3, 0,0,0,0,0,0)
title2use = "hM4Di"
}
data_H$facet_var = "H"
data_FR$facet_var = "Firing Rate"
} else{
if (exp=="camk"){
data_H$score = c(0,0,0,3,3,3, 0,0,0,0,0,0)
data_FR$score = c(0,0,0,3,3,3, 0,0,0,0,0,0)
title2use = "CamkII-hM3D(Gq)"
}else if (exp=="hm4di"){
data_H$score = c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
data_FR$score = c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
title2use = "hM4Di"
}
data_H$facet_var = dtype
data_FR$facet_var = "Firing Rate"
}
data = rbind(data_FR,data_H)
data$facet_var = factor(data$facet_var, levels = c("Firing Rate", dtype))
data$grp = factor(data$grp, levels = c("DREADD","SHAM"))
data$condition = factor(data$condition, levels = c("baseline","transition","treatment"))
p = ggplot(data = data, aes(x = time, y = score, colour = condition)) + facet_grid(. ~ facet_var) +
geom_hline(yintercept = 0) +
geom_line(aes(linetype=grp),size = lineSize) +
scale_linetype_manual(values = c("solid", "dashed")) +
# xlim(xlimits[1],xlimits[2]) +
scale_y_continuous(limits = ylimits) +
xlab("Time Window") + ylab("Baseline Normalized Value") +
geom_vline(xintercept = 15) +
geom_vline(xintercept = 35) +
guides(colour="none", linetype="none") +
ggtitle(title2use) + easy_center_title() +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, size=fontSize),
strip.text.y = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize-2),
axis.title.x = element_text(size=fontSize),
axis.title.y = element_text(size=fontSize),
strip.text = element_text(size = fontSize))
return(p)
}
# run_lme
run_lme <- function(data, dv_var){
# formula
form2use = as.formula(sprintf("%s ~ grp*condition + (1|id)", dv_var))
# run the model
mod2use = lmer(formula = form2use, data = data)
aov_tab = anova(mod2use)
pairwise_res = emmeans(mod2use, pairwise ~ grp*condition)
result = list(model = mod2use,
aov_tab = aov_tab,
pairwise_res = pairwise_res)
return(result)
} # function run_lme
Theoretical plots
CamkII-hM3D(Gq)
camk_plot = make_blankplot(exp = "camk", dtype = "H", fontSize=24)
ggsave(filename = file.path(plotdir,"camk_H_theoretical_prediction.pdf"), plot = camk_plot)
camk_plot
camk_plot = make_blankplot(exp = "camk", dtype = "1/f slope", fontSize=24)
ggsave(filename = file.path(plotdir,"camk_slope_theoretical_prediction.pdf"), plot = camk_plot)
camk_plot
camk_plot = make_blankplot(exp = "camk", dtype = "Broadband Power", fontSize=24)
ggsave(filename = file.path(plotdir,"camk_power_theoretical_prediction.pdf"), plot = camk_plot)
camk_plot
hM4Di
hm4di_plot = make_blankplot(exp = "hm4di", dtype = "H", fontSize=24)
ggsave(filename = file.path(plotdir,"hm4di_H_theoretical_prediction.pdf"), plot = hm4di_plot)
hm4di_plot
hm4di_plot = make_blankplot(exp = "hm4di", dtype = "1/f slope", fontSize=24)
ggsave(filename = file.path(plotdir,"hm4di_slope_theoretical_prediction.pdf"), plot = hm4di_plot)
hm4di_plot
hm4di_plot = make_blankplot(exp = "hm4di", dtype = "Broadband Power", fontSize=24)
ggsave(filename = file.path(plotdir,"hm4di_power_theoretical_prediction.pdf"), plot = hm4di_plot)
hm4di_plot
CamkII-hM3D(Gq) Firing Rate
data_file = here("experimental_data","camk_databasenorm.csv")
data = read.csv(data_file)
data$grp[data$grp=="exp"] = "DREADD"
data$grp[data$grp=="ctrl"] = "SHAM"
data$grp = factor(data$grp)
data2use = data %>%
filter(condition=="treatment")
# Firing Rate
p0 = make_timeplot(data4plot=data, x_var="time", y_var="r_basenorm", grp_var="id",
yLabel="Firing Rate",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_FR_timeplot.pdf"), plot = p0)
p0
# run LME
res = run_lme(data = data, dv_var = "r_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 21.56 21.559 1 11.02 7.3291 0.02036 *
## condition 247.81 123.906 2 763.00 42.1230 < 2e-16 ***
## grp:condition 378.05 189.026 2 763.00 64.2610 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 4.440802e-18
res$aov_tab[3,6]
## [1] 1.613408e-26
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.000 0.644 12.7 -1.39 1.393
## SHAM baseline 0.000 0.509 12.7 -1.10 1.102
## DREADD transition 2.759 0.636 12.2 1.38 4.143
## SHAM transition -0.147 0.503 12.2 -1.24 0.947
## DREADD treatment 3.122 0.631 11.8 1.74 4.500
## SHAM treatment -0.399 0.499 11.8 -1.49 0.690
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.000 0.821 12.7 0.000 1.0000
## DREADD baseline - DREADD transition -2.759 0.262 763.0 -10.531 <.0001
## DREADD baseline - SHAM transition 0.147 0.817 12.5 0.180 1.0000
## DREADD baseline - DREADD treatment -3.122 0.251 763.0 -12.462 <.0001
## DREADD baseline - SHAM treatment 0.399 0.815 12.4 0.490 0.9957
## SHAM baseline - DREADD transition -2.759 0.815 12.4 -3.387 0.0463
## SHAM baseline - SHAM transition 0.147 0.207 763.0 0.710 0.9808
## SHAM baseline - DREADD treatment -3.122 0.811 12.2 -3.850 0.0215
## SHAM baseline - SHAM treatment 0.399 0.198 763.0 2.016 0.3339
## DREADD transition - SHAM transition 2.906 0.811 12.2 3.584 0.0338
## DREADD transition - DREADD treatment -0.363 0.230 763.0 -1.577 0.6142
## DREADD transition - SHAM treatment 3.158 0.809 12.0 3.906 0.0198
## SHAM transition - DREADD treatment -3.269 0.807 11.9 -4.050 0.0157
## SHAM transition - SHAM treatment 0.252 0.182 763.0 1.387 0.7349
## DREADD treatment - SHAM treatment 3.521 0.805 11.8 4.375 0.0093
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
CamkII-hM3D(Gq) Hurst exponent
# H
p1 = make_timeplot(data4plot=data, x_var="time", y_var="H_basenorm", grp_var="id",
yLabel="H",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_H_timeplot.pdf"), plot = p1)
p1
# run LME
res = run_lme(data = data, dv_var = "H_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 13.02 13.020 1 11.05 7.1671 0.02145 *
## condition 116.98 58.489 2 763.00 32.1966 3.767e-14 ***
## grp:condition 195.11 97.556 2 763.00 53.7024 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 3.766969e-14
res$aov_tab[3,6]
## [1] 1.511591e-22
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.000 0.383 14.3 -0.820 0.820
## SHAM baseline 0.000 0.303 14.3 -0.649 0.649
## DREADD transition -0.694 0.375 13.2 -1.503 0.116
## SHAM transition 0.441 0.297 13.2 -0.200 1.081
## DREADD treatment -2.145 0.370 12.5 -2.949 -1.342
## SHAM treatment 0.398 0.293 12.5 -0.237 1.033
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.0000 0.489 14.3 0.000 1.0000
## DREADD baseline - DREADD transition 0.6936 0.206 763.0 3.369 0.0103
## DREADD baseline - SHAM transition -0.4405 0.485 13.9 -0.909 0.9378
## DREADD baseline - DREADD treatment 2.1454 0.197 763.0 10.898 <.0001
## DREADD baseline - SHAM treatment -0.3980 0.482 13.6 -0.825 0.9577
## SHAM baseline - DREADD transition 0.6936 0.482 13.6 1.438 0.7055
## SHAM baseline - SHAM transition -0.4405 0.163 763.0 -2.707 0.0751
## SHAM baseline - DREADD treatment 2.1454 0.479 13.2 4.483 0.0062
## SHAM baseline - SHAM treatment -0.3980 0.156 763.0 -2.557 0.1094
## DREADD transition - SHAM transition -1.1341 0.478 13.2 -2.370 0.2344
## DREADD transition - DREADD treatment 1.4519 0.181 763.0 8.029 <.0001
## DREADD transition - SHAM treatment -1.0915 0.476 12.9 -2.293 0.2638
## SHAM transition - DREADD treatment 2.5860 0.475 12.8 5.448 0.0013
## SHAM transition - SHAM treatment 0.0426 0.143 763.0 0.298 0.9997
## DREADD treatment - SHAM treatment -2.5434 0.472 12.5 -5.386 0.0015
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
p2 = make_scatterplot(data4plot = data2use,
x_var = "r_basenorm",
y_var = "H_basenorm",
grp_var = "grp",
cols2use = cols2use,
xLabel = "Firing Rate",
yLabel = "H",
title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_H_scatterplot.pdf"), plot = p2)
p2
res = cor.test(data2use$r_basenorm, data2use$H_basenorm)
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm and data2use$H_basenorm
## t = -10.554, df = 323, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5830460 -0.4207607
## sample estimates:
## cor
## -0.5063734
res$p.value
## [1] 1.441197e-22
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"],
data2use$H_basenorm[data2use$grp=="DREADD"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$H_basenorm[data2use$grp == "DREADD"]
## t = -3.1517, df = 123, p-value = 0.002039
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4284003 -0.1026790
## sample estimates:
## cor
## -0.2733575
res$p.value
## [1] 0.002038969
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"],
data2use$H_basenorm[data2use$grp=="SHAM"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$H_basenorm[data2use$grp == "SHAM"]
## t = -3.3387, df = 198, p-value = 0.001005
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.35813380 -0.09517126
## sample estimates:
## cor
## -0.2308639
res$p.value
## [1] 0.001005481
CamkII-hM3D(Gq) 1/f slope
# 1/f slope
p3 = make_timeplot(data4plot=data, x_var="time", y_var="slope_basenorm", grp_var="id",
yLabel="1/f slope",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_slope_timeplot.pdf"), plot = p3)
p3
# run LME
res = run_lme(data = data, dv_var = "slope_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 0.6579 0.65793 1 11.11 0.6340 0.4426
## condition 0.6974 0.34872 2 763.00 0.3360 0.7147
## grp:condition 4.0365 2.01827 2 763.00 1.9447 0.1437
res$aov_tab[2,6]
## [1] 0.7147202
res$aov_tab[3,6]
## [1] 0.1437322
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.0000 0.201 19.9 -0.419 0.419
## SHAM baseline 0.0000 0.159 19.9 -0.331 0.331
## DREADD transition 0.1994 0.192 16.7 -0.206 0.605
## SHAM transition -0.1811 0.152 16.7 -0.502 0.140
## DREADD treatment 0.0171 0.186 14.8 -0.381 0.415
## SHAM treatment -0.1305 0.147 14.8 -0.445 0.184
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.0000 0.256 19.9 0.000 1.0000
## DREADD baseline - DREADD transition -0.1994 0.156 763.0 -1.281 0.7954
## DREADD baseline - SHAM transition 0.1811 0.252 18.6 0.720 0.9770
## DREADD baseline - DREADD treatment -0.0171 0.149 763.0 -0.115 1.0000
## DREADD baseline - SHAM treatment 0.1305 0.249 17.9 0.524 0.9944
## SHAM baseline - DREADD transition -0.1994 0.249 17.9 -0.800 0.9637
## SHAM baseline - SHAM transition 0.1811 0.123 763.0 1.472 0.6822
## SHAM baseline - DREADD treatment -0.0171 0.245 16.7 -0.070 1.0000
## SHAM baseline - SHAM treatment 0.1305 0.118 763.0 1.110 0.8774
## DREADD transition - SHAM transition 0.3805 0.245 16.7 1.555 0.6364
## DREADD transition - DREADD treatment 0.1823 0.137 763.0 1.334 0.7662
## DREADD transition - SHAM treatment 0.3299 0.242 15.9 1.363 0.7471
## SHAM transition - DREADD treatment -0.1982 0.240 15.5 -0.824 0.9584
## SHAM transition - SHAM treatment -0.0506 0.108 763.0 -0.468 0.9972
## DREADD treatment - SHAM treatment 0.1476 0.238 14.8 0.621 0.9876
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
p4 = make_scatterplot(data4plot = data2use,
x_var = "r_basenorm",
y_var = "slope_basenorm",
grp_var = "grp",
cols2use = cols2use,
xLabel = "Firing Rate",
yLabel = "1/f slope",
title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_slope_scatterplot.pdf"), plot = p4)
p4
res = cor.test(data2use$r_basenorm, data2use$slope_basenorm)
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm and data2use$slope_basenorm
## t = 1.3504, df = 323, p-value = 0.1778
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.03414132 0.18223593
## sample estimates:
## cor
## 0.07492923
res$p.value
## [1] 0.1778205
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"],
data2use$slope_basenorm[data2use$grp=="DREADD"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$slope_basenorm[data2use$grp == "DREADD"]
## t = 2.3733, df = 123, p-value = 0.01918
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.03492914 0.37121966
## sample estimates:
## cor
## 0.2092531
res$p.value
## [1] 0.01918032
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"],
data2use$slope_basenorm[data2use$grp=="SHAM"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$slope_basenorm[data2use$grp == "SHAM"]
## t = -0.078323, df = 198, p-value = 0.9377
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1441957 0.1332779
## sample estimates:
## cor
## -0.005566076
res$p.value
## [1] 0.9376503
CamkII-hM3D(Gq) Broadband power
# broadband power
p5 = make_timeplot(data4plot=data, x_var="time", y_var="isp_basenorm", grp_var="id",
yLabel="Power",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_power_timeplot.pdf"), plot = p5)
p5
# run LME
res = run_lme(data = data, dv_var = "isp_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 4.529 4.529 1 11.1 3.5611 0.08557 .
## condition 121.880 60.940 2 763.0 47.9183 < 2.2e-16 ***
## grp:condition 37.754 18.877 2 763.0 14.8435 4.74e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 2.490172e-20
res$aov_tab[3,6]
## [1] 4.740382e-07
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.0000 0.236 18.5 -0.494 0.4940
## SHAM baseline 0.0000 0.186 18.5 -0.391 0.3906
## DREADD transition -0.4196 0.226 15.8 -0.900 0.0609
## SHAM transition -0.0173 0.179 15.8 -0.397 0.3626
## DREADD treatment -1.4785 0.221 14.3 -1.951 -1.0059
## SHAM treatment -0.3848 0.175 14.3 -0.758 -0.0112
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.0000 0.300 18.5 0.000 1.0000
## DREADD baseline - DREADD transition 0.4196 0.172 763.0 2.436 0.1453
## DREADD baseline - SHAM transition 0.0173 0.296 17.4 0.058 1.0000
## DREADD baseline - DREADD treatment 1.4785 0.165 763.0 8.976 <.0001
## DREADD baseline - SHAM treatment 0.3848 0.293 16.8 1.313 0.7745
## SHAM baseline - DREADD transition 0.4196 0.293 16.8 1.431 0.7089
## SHAM baseline - SHAM transition 0.0173 0.136 763.0 0.127 1.0000
## SHAM baseline - DREADD treatment 1.4785 0.289 15.8 5.119 0.0012
## SHAM baseline - SHAM treatment 0.3848 0.130 763.0 2.955 0.0378
## DREADD transition - SHAM transition -0.4023 0.289 15.8 -1.394 0.7301
## DREADD transition - DREADD treatment 1.0589 0.151 763.0 6.999 <.0001
## DREADD transition - SHAM treatment -0.0348 0.286 15.2 -0.122 1.0000
## SHAM transition - DREADD treatment 1.4612 0.284 14.8 5.142 0.0014
## SHAM transition - SHAM treatment 0.3675 0.120 763.0 3.072 0.0266
## DREADD treatment - SHAM treatment -1.0937 0.281 14.3 -3.887 0.0161
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
p6 = make_scatterplot(data4plot = data2use,
x_var = "r_basenorm",
y_var = "isp_basenorm",
grp_var = "grp",
cols2use = cols2use,
xLabel = "Firing Rate",
yLabel = "Power",
title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_power_scatterplot.pdf"), plot = p6)
p6
res = cor.test(data2use$r_basenorm, data2use$isp_basenorm)
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm and data2use$isp_basenorm
## t = -3.8348, df = 323, p-value = 0.0001511
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3104191 -0.1022020
## sample estimates:
## cor
## -0.2086741
res$p.value
## [1] 0.0001511161
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"],
data2use$isp_basenorm[data2use$grp=="DREADD"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$isp_basenorm[data2use$grp == "DREADD"]
## t = -0.6226, df = 123, p-value = 0.5347
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2293995 0.1207461
## sample estimates:
## cor
## -0.05604991
res$p.value
## [1] 0.5346984
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"],
data2use$isp_basenorm[data2use$grp=="SHAM"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$isp_basenorm[data2use$grp == "SHAM"]
## t = 1.6355, df = 198, p-value = 0.1035
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.0236671 0.2501865
## sample estimates:
## cor
## 0.115453
res$p.value
## [1] 0.103532
hM4Di Firing Rate
data_file = here("experimental_data","hm4d_databasenorm.csv")
data = read.csv(data_file)
data$grp[data$grp=="exp"] = "DREADD"
data$grp[data$grp=="ctrl"] = "SHAM"
data$grp = factor(data$grp)
data2use = data %>%
filter(condition=="treatment")
# Firing Rate
p0 = make_timeplot(data4plot=data, x_var="time", y_var="r_basenorm", grp_var="id",
yLabel="Firing Rate",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_FR_timeplot.pdf"), plot = p0)
p0
# run LME
res = run_lme(data = data, dv_var = "r_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 4.959 4.959 1 8.01 2.5621 0.1481
## condition 157.418 78.709 2 586.00 40.6677 <2e-16 ***
## grp:condition 244.156 122.078 2 586.00 63.0757 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 2.891814e-17
res$aov_tab[3,6]
## [1] 1.54989e-25
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.000 0.692 8.69 -1.57 1.573
## SHAM baseline 0.000 0.692 8.69 -1.57 1.573
## DREADD transition -1.777 0.687 8.46 -3.35 -0.208
## SHAM transition -0.337 0.687 8.46 -1.91 1.232
## DREADD treatment -2.836 0.684 8.32 -4.40 -1.268
## SHAM treatment 0.328 0.684 8.32 -1.24 1.895
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.000 0.978 8.69 0.000 1.0000
## DREADD baseline - DREADD transition 1.777 0.213 586.00 8.363 <.0001
## DREADD baseline - SHAM transition 0.337 0.975 8.57 0.346 0.9991
## DREADD baseline - DREADD treatment 2.836 0.203 586.00 13.956 <.0001
## DREADD baseline - SHAM treatment -0.328 0.973 8.50 -0.337 0.9992
## SHAM baseline - DREADD transition 1.777 0.975 8.57 1.823 0.4994
## SHAM baseline - SHAM transition 0.337 0.213 586.00 1.587 0.6073
## SHAM baseline - DREADD treatment 2.836 0.973 8.50 2.914 0.1286
## SHAM baseline - SHAM treatment -0.328 0.203 586.00 -1.612 0.5908
## DREADD transition - SHAM transition -1.440 0.972 8.46 -1.482 0.6833
## DREADD transition - DREADD treatment 1.059 0.187 586.00 5.671 <.0001
## DREADD transition - SHAM treatment -2.105 0.970 8.39 -2.171 0.3394
## SHAM transition - DREADD treatment 2.498 0.970 8.39 2.577 0.2037
## SHAM transition - SHAM treatment -0.665 0.187 586.00 -3.562 0.0053
## DREADD treatment - SHAM treatment -3.163 0.968 8.32 -3.269 0.0805
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
hM4Di Hurst exponent
# H
p1 = make_timeplot(data4plot=data, x_var="time", y_var="H_basenorm", grp_var="id",
yLabel="H",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_H_timeplot.pdf"), plot = p1)
p1
# run LME
res = run_lme(data = data, dv_var = "H_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 4.36 4.362 1 8.01 1.779 0.219
## condition 320.59 160.294 2 586.00 65.376 <2e-16 ***
## grp:condition 196.79 98.394 2 586.00 40.130 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 2.348727e-26
res$aov_tab[3,6]
## [1] 4.637511e-17
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.000 0.867 8.55 -1.977 1.98
## SHAM baseline 0.000 0.867 8.55 -1.977 1.98
## DREADD transition 2.089 0.862 8.36 0.116 4.06
## SHAM transition 0.157 0.862 8.36 -1.816 2.13
## DREADD treatment 3.296 0.859 8.25 1.325 5.27
## SHAM treatment 0.402 0.859 8.25 -1.569 2.37
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.000 1.226 8.55 0.000 1.0000
## DREADD baseline - DREADD transition -2.089 0.239 586.00 -8.733 <.0001
## DREADD baseline - SHAM transition -0.157 1.223 8.46 -0.128 1.0000
## DREADD baseline - DREADD treatment -3.296 0.229 586.00 -14.410 <.0001
## DREADD baseline - SHAM treatment -0.402 1.221 8.40 -0.329 0.9993
## SHAM baseline - DREADD transition -2.089 1.223 8.46 -1.708 0.5601
## SHAM baseline - SHAM transition -0.157 0.239 586.00 -0.655 0.9866
## SHAM baseline - DREADD treatment -3.296 1.221 8.40 -2.700 0.1729
## SHAM baseline - SHAM treatment -0.402 0.229 586.00 -1.757 0.4945
## DREADD transition - SHAM transition 1.932 1.219 8.36 1.585 0.6275
## DREADD transition - DREADD treatment -1.207 0.210 586.00 -5.745 <.0001
## DREADD transition - SHAM treatment 1.687 1.217 8.31 1.386 0.7346
## SHAM transition - DREADD treatment -3.139 1.217 8.31 -2.579 0.2039
## SHAM transition - SHAM treatment -0.245 0.210 586.00 -1.167 0.8524
## DREADD treatment - SHAM treatment 2.894 1.215 8.25 2.381 0.2632
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
p2 = make_scatterplot(data4plot = data2use,
x_var = "r_basenorm",
y_var = "H_basenorm",
grp_var = "grp",
cols2use = cols2use,
xLabel = "Firing Rate",
yLabel = "H",
title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_H_scatterplot.pdf"), plot = p2)
p2
res = cor.test(data2use$r_basenorm, data2use$H_basenorm)
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm and data2use$H_basenorm
## t = -13.91, df = 248, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7264211 -0.5860897
## sample estimates:
## cor
## -0.6620185
res$p.value
## [1] 6.649365e-33
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"],
data2use$H_basenorm[data2use$grp=="DREADD"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$H_basenorm[data2use$grp == "DREADD"]
## t = -14.142, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8456471 -0.7093022
## sample estimates:
## cor
## -0.786895
res$p.value
## [1] 1.481409e-27
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"],
data2use$H_basenorm[data2use$grp=="SHAM"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$H_basenorm[data2use$grp == "SHAM"]
## t = -9.99, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7560316 -0.5594230
## sample estimates:
## cor
## -0.6692811
res$p.value
## [1] 1.440809e-17
hM4Di 1/f slope
# 1/f slope
p3 = make_timeplot(data4plot=data, x_var="time", y_var="slope_basenorm", grp_var="id",
yLabel="1/f slope",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_slope_timeplot.pdf"), plot = p3)
p3
# run LME
res = run_lme(data = data, dv_var = "slope_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 3.9455 3.9455 1 8.06 2.9629 0.1232
## condition 0.3827 0.1913 2 586.00 0.1437 0.8662
## grp:condition 28.3783 14.1891 2 586.00 10.6556 2.848e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 0.8661809
res$aov_tab[3,6]
## [1] 2.847621e-05
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.000 0.257 12.54 -0.5575 0.557
## SHAM baseline 0.000 0.257 12.54 -0.5575 0.557
## DREADD transition -0.358 0.248 10.91 -0.9047 0.189
## SHAM transition 0.229 0.248 10.91 -0.3181 0.776
## DREADD treatment -0.570 0.243 9.99 -1.1114 -0.029
## SHAM treatment 0.524 0.243 9.99 -0.0171 1.065
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.000 0.364 12.54 0.000 1.0000
## DREADD baseline - DREADD transition 0.358 0.176 586.00 2.029 0.3271
## DREADD baseline - SHAM transition -0.229 0.357 11.71 -0.640 0.9853
## DREADD baseline - DREADD treatment 0.570 0.169 586.00 3.383 0.0099
## DREADD baseline - SHAM treatment -0.524 0.354 11.23 -1.482 0.6814
## SHAM baseline - DREADD transition 0.358 0.357 11.71 1.001 0.9088
## SHAM baseline - SHAM transition -0.229 0.176 586.00 -1.299 0.7859
## SHAM baseline - DREADD treatment 0.570 0.354 11.23 1.612 0.6074
## SHAM baseline - SHAM treatment -0.524 0.169 586.00 -3.109 0.0240
## DREADD transition - SHAM transition -0.587 0.351 10.91 -1.670 0.5749
## DREADD transition - DREADD treatment 0.213 0.155 586.00 1.373 0.7434
## DREADD transition - SHAM treatment -0.882 0.347 10.45 -2.539 0.1969
## SHAM transition - DREADD treatment 0.799 0.347 10.45 2.301 0.2740
## SHAM transition - SHAM treatment -0.295 0.155 586.00 -1.906 0.3991
## DREADD treatment - SHAM treatment -1.094 0.343 9.99 -3.186 0.0774
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
p4 = make_scatterplot(data4plot = data2use,
x_var = "r_basenorm",
y_var = "slope_basenorm",
grp_var = "grp",
cols2use = cols2use,
xLabel = "Firing Rate",
yLabel = "1/f slope",
title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_slope_scatterplot.pdf"), plot = p4)
p4
res = cor.test(data2use$r_basenorm, data2use$slope_basenorm)
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm and data2use$slope_basenorm
## t = 6.6842, df = 248, p-value = 1.523e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2802265 0.4909770
## sample estimates:
## cor
## 0.3907097
res$p.value
## [1] 1.523323e-10
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"],
data2use$slope_basenorm[data2use$grp=="DREADD"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$slope_basenorm[data2use$grp == "DREADD"]
## t = 4.898, df = 123, p-value = 2.98e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2458295 0.5412082
## sample estimates:
## cor
## 0.4039967
res$p.value
## [1] 2.980145e-06
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"],
data2use$slope_basenorm[data2use$grp=="SHAM"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$slope_basenorm[data2use$grp == "SHAM"]
## t = 2.1629, df = 123, p-value = 0.03249
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.01635512 0.35508475
## sample estimates:
## cor
## 0.1914129
res$p.value
## [1] 0.03248631
hM4Di Broadband power
# broadband power
p5 = make_timeplot(data4plot=data, x_var="time", y_var="isp_basenorm", grp_var="id",
yLabel="Power",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_power_timeplot.pdf"), plot = p5)
p5
# run LME
res = run_lme(data = data, dv_var = "isp_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## grp 25.67 25.67 1 8.02 1.4412 0.2642
## condition 2073.29 1036.64 2 586.00 58.1915 < 2.2e-16 ***
## grp:condition 409.89 204.95 2 586.00 11.5045 1.257e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 8.867731e-24
res$aov_tab[3,6]
## [1] 1.256753e-05
res$pairwise_res
## $emmeans
## grp condition emmean SE df lower.CL upper.CL
## DREADD baseline 0.000 1.53 9.36 -3.448 3.45
## SHAM baseline 0.000 1.53 9.36 -3.448 3.45
## DREADD transition 4.000 1.51 8.90 0.570 7.43
## SHAM transition 0.498 1.51 8.90 -2.933 3.93
## DREADD treatment 6.643 1.50 8.62 3.223 10.06
## SHAM treatment 2.634 1.50 8.62 -0.787 6.05
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## DREADD baseline - SHAM baseline 0.000 2.168 9.36 0.000 1.0000
## DREADD baseline - DREADD transition -4.000 0.645 586.00 -6.205 <.0001
## DREADD baseline - SHAM transition -0.498 2.155 9.13 -0.231 0.9999
## DREADD baseline - DREADD treatment -6.643 0.616 586.00 -10.776 <.0001
## DREADD baseline - SHAM treatment -2.634 2.146 8.99 -1.227 0.8143
## SHAM baseline - DREADD transition -4.000 2.155 9.13 -1.857 0.4795
## SHAM baseline - SHAM transition -0.498 0.645 586.00 -0.772 0.9722
## SHAM baseline - DREADD treatment -6.643 2.146 8.99 -3.095 0.0961
## SHAM baseline - SHAM treatment -2.634 0.616 586.00 -4.272 0.0003
## DREADD transition - SHAM transition 3.503 2.141 8.90 1.636 0.5981
## DREADD transition - DREADD treatment -2.643 0.566 586.00 -4.667 0.0001
## DREADD transition - SHAM treatment 1.367 2.132 8.76 0.641 0.9843
## SHAM transition - DREADD treatment -6.145 2.132 8.76 -2.882 0.1320
## SHAM transition - SHAM treatment -2.136 0.566 586.00 -3.772 0.0024
## DREADD treatment - SHAM treatment 4.009 2.124 8.62 1.888 0.4663
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 6 estimates
p6 = make_scatterplot(data4plot = data2use,
x_var = "r_basenorm",
y_var = "isp_basenorm",
grp_var = "grp",
cols2use = cols2use,
xLabel = "Firing Rate",
yLabel = "Power",
title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_power_scatterplot.pdf"), plot = p6)
p6
res = cor.test(data2use$r_basenorm, data2use$isp_basenorm)
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm and data2use$isp_basenorm
## t = -8.017, df = 248, p-value = 4.275e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5469561 -0.3492670
## sample estimates:
## cor
## -0.4536751
res$p.value
## [1] 4.275038e-14
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"],
data2use$isp_basenorm[data2use$grp=="DREADD"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$isp_basenorm[data2use$grp == "DREADD"]
## t = -6.0163, df = 123, p-value = 1.892e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6020279 -0.3287518
## sample estimates:
## cor
## -0.4768312
res$p.value
## [1] 1.892143e-08
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"],
data2use$isp_basenorm[data2use$grp=="SHAM"])
res
##
## Pearson's product-moment correlation
##
## data: data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$isp_basenorm[data2use$grp == "SHAM"]
## t = -4.4494, df = 123, p-value = 1.904e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5143217 -0.2105001
## sample estimates:
## cor
## -0.3723438
res$p.value
## [1] 1.903682e-05
Make figures of power spectrum
freq_data = read.csv(here("experimental_data","freq.csv"))
freq_names = freq_data$colnames2use
# CamkII=-hM3D(Gq)
data_file = here("experimental_data","LFPspectra_xsub_camk.csv")
data = read.csv(data_file)
mean_data = data %>% filter(condition=="treatment")
# n1 = length(unique(mean_data$id[mean_data$grp=="SHAM"]))
n1 = sum(mean_data$grp=="SHAM")
data1 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="SHAM",freq_names]), grp = "SHAM", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="SHAM",freq_names], 2, sd) / sqrt(n1)))
# n2 = length(unique(mean_data$id[mean_data$grp=="DREADD"]))
n2 = sum(mean_data$grp=="DREADD")
data2 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="DREADD",freq_names]), grp = "DREADD", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="DREADD",freq_names], 2, sd) / sqrt(n2)))
data2plot = rbind(data1,data2)
data2plot$logFreq = log10(data2plot$freq)
camk_pow_spec = ggplot(data = data2plot,
aes(x = logFreq,
y = mean_spectrum_bn,
colour=grp,
fill=grp)) +
geom_ribbon(aes(ymin = mean_spectrum_bn - se_spectrum_bn,
ymax = mean_spectrum_bn + se_spectrum_bn),
alpha = 0.2) +
geom_line() +
geom_hline(yintercept=0) +
xlab("Frequency (Hz)") +
ylab("Baseline Normalized Power") +
scale_x_continuous(labels = c(0.1, 1, 10, 100)) +
ggtitle("CamkII-hM3D(Gq)") + easy_center_title() +
easy_add_legend_title("Group") + easy_adjust_legend(to = "center") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, size=fontSize),
strip.text.y = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize-2),
axis.title.x = element_text(size=fontSize),
axis.title.y = element_text(size=fontSize),
strip.text = element_text(size = fontSize))
# hM4Di
data_file = here("experimental_data","LFPspectra_xsub_hm4d.csv")
data = read.csv(data_file)
mean_data = data %>% filter(condition=="treatment")
# n1 = length(unique(mean_data$id[mean_data$grp=="SHAM"]))
n1 = sum(mean_data$grp=="SHAM")
data1 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="SHAM",freq_names]), grp = "SHAM", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="SHAM",freq_names], 2, sd) / sqrt(n1)))
# n2 = length(unique(mean_data$id[mean_data$grp=="DREADD"]))
n2 = sum(mean_data$grp=="DREADD")
data2 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="DREADD",freq_names]), grp = "DREADD", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="DREADD",freq_names], 2, sd) / sqrt(n2)))
data2plot = rbind(data1,data2)
data2plot$logFreq = log10(data2plot$freq)
hm4di_pow_spec = ggplot(data = data2plot,
aes(x = logFreq,
y = mean_spectrum_bn,
colour=grp,
fill=grp)) +
geom_ribbon(aes(ymin = mean_spectrum_bn - se_spectrum_bn,
ymax = mean_spectrum_bn + se_spectrum_bn),
alpha = 0.2) +
geom_line() +
geom_hline(yintercept=0) +
xlab("Frequency (Hz)") +
ylab("Baseline Normalized Power") +
scale_x_continuous(labels = c(0.1, 1, 10, 100)) +
ggtitle("hM4Di") + easy_center_title() +
easy_add_legend_title("Group") + easy_adjust_legend(to = "center") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, size=fontSize),
strip.text.y = element_text(size=fontSize),
axis.text.x = element_text(size=fontSize),
axis.text.y = element_text(size=fontSize-2),
axis.title.x = element_text(size=fontSize),
axis.title.y = element_text(size=fontSize),
strip.text = element_text(size = fontSize))
p_final = camk_pow_spec / hm4di_pow_spec + plot_layout(guides = "collect")
ggsave(filename = file.path(plotdir, "dreadd_pow_spec.pdf"),
plot = p_final,
width = 8,height=8)
p_final
